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1 Introduction 



One frequently comes across complex valued path integral weights in physics. 
A major example is the study of field theories in Minkowski space. Even in 
^ ' Euclidean field theories the effective action may turn out to be complex, e.g. 

■ QCD with a chemical potential. Another case, in which we will be interested 

in this paper, is the study of the solution space of Schwinger-Dyson equations 
as the phase space of the associated quantum field theory [II2]- 

Standard Monte Carlo methods do not work in the case, where the path 
integral weight is not positive-definite. This is the famous "sign problem". 
The root of the problem lies in the fact that Monte Carlo methods work on 
the probabilistic interpretation of the path integral weight , where S is the 
action of the system. In the examples above, the exponent becomes complex 
and the probabilistic interpretation fails. 
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Complex Langevin equations were proposed by Parisi [3] and Klauder [1] to 
simulate systems with complex valued path integral weights. For d dimensional 
systems with real weights, a Langevin equation m d + 1 dimensions may be 
used to study the partition function of the system. When properly set up, the 
Langevin process converges to a unique stationary distribution, which is the 
partition function of the associated system, in the limit of the additional di- 
mension going to infinity. This fact was used by Parisi and Wu in the stochastic 
quantization of quantum fields |5j . For systems with complex actions, one can 
still write down a (complex) Langevin equation, as suggested by Parisi and 
Klauder, but this approach comes with many problems. First of all, it is not 
certain that the complex Langevin simulation will ever converge to a stationary 
distribution and if it does, there may be many such stationary distributions, 
see e.g. j6H71l8|9l[T0] . Salcedo noted that these stationary distributions may be 
constructed by path integrals over contours that connect zeros of the path in- 
tegral weight e""^ [H] . Other authors noted that these stationary distributions 
satisfy Schwinger-Dyson equations, e.g. [Tl|12|13pl] . On a completely differ- 
ent track of research, one of the authors with Garcia and Guralnik studied the 
solution space of Schwinger-Dyson equations and noted that different solutions 
to Schwinger-Dyson equations may be written as path integrals over contours 
that connect zeros of the path integral weight [US], exactly as Salcedo 
suggested for stationary distributions of the complex Langevin equation. Fur- 
thermore, the authors of references [Ij and [2] suggested different solutions of 
the Schwinger-Dyson equations be interpreted as different phases of the as- 
sociated quantum field theory and studied this proposal in detail. Salcedo in 
his paper [H] made the same suggestion about stationary distributions of the 
complex Langevin equations in one sentence at the conclusion, but did not 
elaborate on it. Our aim in this paper is to point out and clarify the connec- 
tion between these two lines of research and propose @@ complex Langevin 
equations as a basis of a new numerical method of studying different phases 
of a quantum field theory. We start by a rather detailed explanation of the 
mechanism of complex Langevin equations and the problems associated with 
them in section [2l Section [3] contains the main result of this paper. There, in a 
zero dimensional setting, we show that the stationary distributions of a com- 
plex Langevin equation are the solutions of the Schwinger-Dyson equations 
for the associated quantum field theory. Furthermore, these solutions may be 
constructed by changing the integration contour of path integrals from real 
paths to contours that connect the zeros of the path integral weight e""^ on the 
complex plane. In section HI we do the trivial generalization of the problem to 
a lattice and discuss related issues. We conclude by further summarizing the 
results of [T] and [2] and point out the connection between different phases 
of a quantum field theory and stationary distributions of complex Langevin 
equations. Based on this observation, we propose complex Langevin equations 
as a numerical method of studying different phases of a quantum field theory. 
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2 The Complex Langevin Equation 



In this section, we introduce Langevin equations. We discuss zero dimensional 
quantum field theories with a scalar field for simplicity. Generalization to 
vector fields and higher dimensions is straightforward, see e.g. [15J. For systems 
with action S {(f)), we are interested in calculating expectation values like 

using Langevin equations. We will assume 5* (z) to be an analytic function of 
the complex variable z. 

Note that our formulation is a Euclidean space formulation, i.e. the path in- 
tegral weight will be complex only when the action is complex. However in 
Minkowski space formulations of quantum field theories, the path integral 
weight will still be complex even though the action is real, since the weight 
is defined by e'^^'^^^'^\ When we speak about complex actions in Euclidean 
space, our results will be applicable to Minkowski space actions after nec- 
essary modifications. One has to note that S{(j)) = —iSM^cj)) and introduce 
appropriate factors of z's in the generating function definition and terms of 
Schwinger-Dyson equation. 



2. 1 Real Actions 



When the action is real (for the moment {/) is a real field), one can create a 
stochastic process using a Langevin equation with a unique stationary distri- 
bution e-^('^V/#e-^('^): 

dS 

d(t){T) = -——dr + dw{r), (2) 
90(r) 

where r is a fictitious time and w (r) is the Wiener process normalized to 
satisfy: 

{dw{T)) = 0, {dw{T)dw{T)) = 2dT, {dw{Ti)dw{T2)) =0 (n ^ rs). 

(3) 

Then one can run this Langevin process to calculate the intended expectation 
values as in equation ([1]). We first show that this Langevin process really 
converges to the intended stationary distribution. 
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Associated with the Langevin process is a probabihty density r) 

(F(0(r)))= /f/0F(0)P(0,r), (4) 



(5) 



which can be shown (e.g. [12]) to satisfy the Fokker-Plack equation 
dPi<l>, r) d ( d dS\ 



Note that the Fokker-Planck equation enables the normahzation condition 

d(j)P{(p,T) = 1 (6) 



to be independent of time since 

d_ 

dr 



d 



d<pP{(t),T) = 

with an appropriate boundary condition on P{(j), r 



(7) 



Now we look at the asymptotic behavior of P{4>, r) as fictitous time goes to 
infinity. Introducing the quantity 



g(0,r) = P(0,r)e^W/^ 
we can rewrite equation ([5]) as 



(8) 



dr 






-H 



FP 



T . 



(9) 



where the Hermitian operator Hpp is called the Fokker-Planck Hamiltonian. 
Then 

Q,(^) = (10) 

is an eigenstate of the Hamiltonian Hpp with eigenvalue 0. Furthermore, (as- 
suming Qq G L^) it is the ground state since it is nowhere vanishing. Then 
the time independent eigenvalue equation 



HFpQn{<P) = EnQn{<P), 



has solutions with the property > 0, for n > 0. Using these eigenvalues 
and eigenf unctions, one can write any solution to equation (Q as 



Q(0, r) = P(0, r)e^(^)/2 = E Qni^^ 



:i2) 



n=0 
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Since we are looking for solutions P{(f),T) that are probability distributions, 
we apply the normalization condition ([6]), which sets Oq = 1/ / d(j)e~^'^'^\ This 
result is obtained by using the orthonormality of the eigenf unctions Qn- Other 
coefficients are set by the initial probability distribution associated with the 
random variable 0(r) at r = 0. Then the limit 

lim Q{(f),T) = aoQo{4>), (13) 

T — 5-00 

implies that one recovers the desired stationary probability distribution for 
the Langevin process as the fictitous time goes to infinity: 

lim P(0,r) = (14) 

Note that this result is independent of the initial conditions. Going back to 
the expectation value problem ([1]), 

Jim (O(0(r))) = (O(0)), (15) 

with {0{(f))) given by ([T]) and ergodicity assures the averaging over the path 
prescription 

(O(0))= hm i r 0{<P{T))dT. (16) 

These observations have been used by Parisi and Wu in the past to formu- 
late the stochastic quantization of quantum fields |5]. See e.g. p!7|15|18] for a 
detailed discussion. 



2.2 Complex Actions 



Now we turn to the case where the action is complex. We want to know if 
we can still use the Langevin equation to calculate desired expectation values. 
We start by rewriting equation ([2]) (now S being complex) in terms of two 
real variables (Pr{t) and 0/(r) as 



;/)/(r) = —Im 



dS 



dS 



dr + dwij-) 



(17) 



where 0(r) = 4>r{t) +i(j)i{T). w{t) is again the Wiener process that is normal- 
ized to satisfy the mean and variance conditions of equation ([3]). Note that the 
equation for d(f)i has a zero diffusion coefficient (see [Ldj for an example where 
it is not zero), but still is a stochastic equation through its dependence on 0^. 
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Complex Langevin equations may be modified to include a term, the kernel, 
that may be useful to stabilize the system, e.g [20f28] . For our purposes, we fo- 
cus on equation fll7l) . Note that we have two different random variables 0_r(t) 
and 0/(t), therefore the real probability distribution associated with equation 
fll7l) will be of the form 



(F(0(r))) = J dc^Rd^j F{(t>R + i<t>i)P{<t>R, h, r) (18) 
where we assume F to be analytic. 

There are two important questions related to this process. The first question 
is whether the probability distribution P(0_r, 0/, r) converges to a stationary 
distribution at all, 

lim P(0^,0,,r) =P(0R,0,). (19) 

T — ^OO 

If it does, does it converge to the desired result: 



d(l)jid(f)i F{(f)R + i(t)i)P{(t>R, (pi 



^Jd^nFiM^-^ 



(20) 



None of these questions have been completely answered so far. Some rigorous 
conditions to verify aposteriori the correct convergence of the process are given 
in 

To understand the difficulties related to the convergence problem, we derive 
the Fokker-Planck equation. First we note that applying the rules of Ito cal- 
culus (see for example [16]) to the complex Langevin equation (fTTj) will give 
the identity 



d_ 



d^F 


dF 


Re 






d<Pl 


d(j)R 




d^F 


OF 


'dS 


502 


50 


50 



dS 



OF 



Im 



dS 



(21) 



where the last line follows from the analyticity of F{(f)). Then using equation 
f|T8l) one can show that P{(f)R, 0/, r) (assuming appropriate differentiability 
and boundary conditions) satisfies the following Fokker-Planck equation: 



5P(0H,0/,r) 
5r 



OfpP{4>r, (Pi, t) 
/ 52 5 „ 



^1 



dS 



dS 



P(0^,0,,r). (22) 



A general statement on the existence of a unique zero eigenvalue (stationary) 
solution P{(j)R,(f)j) for the operator Opp cannot be made. Furthermore, zero 
eigenvalue solutions may exist in the sense of distributions [8] . 
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One can assume a complex valued function P{4>r, t) on the real axis such that 

j^d(l)RF{(t)R)P{(t)R,r) = J d(t)Rd(l)lF{(t)R + l(l)l)P{(t)R,(l)l,T), (23) 

based on the implicit assumption that this equation actually has a solution 
[3]. The reverse question, existence of a positive P{4>r, 4>i, t) given a complex 
P{(f)R,T) is disscussed in [2T|22f23j . Using this definition, analyticity of -F(0) 
and integration by parts in equation (l2Til one can show that P{4>r, t) satisfies 
the pseudo Fokker-Planck equation 

^%ll = 6.pP(fcr) = ^f^ + ^)p(fcr). (24) 

or d(f)R \d(f)R d(f)Rj 

which has the same form as that of equation ([5]). 

A formal solution to equation (l23l) was introduced in [23]. First note that 

F{<f)R + t<f)i) = e^'^F{<f)R), (25) 

due to the analyticity of F{(f)). Inserting this statement into equation (l23ll and 
performing the partial integration assuming necessary boundary and differen- 
tiability conditions, 

P{(I)R, t) = J dcPi e-"^'^P{<j)R, <j)j, t). (26) 



Because the action S{(f)R) is complex, one cannot make general statements 
about the spectrum oi Opp (see [25] for the spectral theorem for a limited 
class of such operators). Furthermore, the relation between the pseudo Fokker- 
Planck equation and the complex Langevin equation were derived based on 
certain assumptions. This should be understood in the sense of distributions, 
meaning only a formal expression of the identity (12T]) . Note that e"'^^'^^^ is 
still a stationary solution (i.e. Opp e~^^^^^ = 0), but in general the stationary 
solutions exist as distributions and the uniqueness of stationary solution is not 
certain [H]. 



3 Stationary Distributions of the Complex Langevin Equation and 
the Boundary Conditions of the Schwinger-Dyson Equation 



Despite the difficulties in proving rigorous results, complex Langevin simula- 
tions have been used to study many different problems. Interesting cases are 
those for which the simulation converges to a stationary distribution which 
is not equivalent to the original complex distribution, e.g. 
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This must be related to the existence of other stationary distributions. Here 
we discuss a conjecture related to the stationary distributions by Salcedo [9] 
and the well know result that stationary distributions satisfy Schwinger-Dyson 
identities, e.g. [TT|12fl3|[T^ . We show that they follow from one another, based 
on other work on the boundary conditions of Schwinger-Dyson equations |T|2j . 
We will point out the results of these references during our discussion. 

We start by assuming that the complex Langevin process has a stationary 
state. For the stationary distribution, the LHS of equation (12T!) must be zero. 
Then setting = 0", n = 1,2,3,..., reproduces the Schwinger-Dyson 

identities for the Green's functions of the quantum field theory defined by the 
action S, i.e. 




n = k ^0-^|^^ = (^-l)(0-^). (27) 

Now define the generating function Z{j) for the stationary distribution 

m = E ^ = {e^') ■ (28) 

We will assume that the radius of convergence for this series is nonzero. Then 
there exists a neighbourhood of j around j = such that the Schwinger-Dyson 
differential equation holds, i.e. 

Zij)=jZij). (29) 

To see that (129|) produces the same identites as (127|) . substitute the definition of 
the generating function Z{j), equation (1251) . in the Schwinger-Dyson equation, 
differentiate with respect to j an appropriate number of times and set j = 
at the end. One gets the identities (l27j) order by order at the end of this 
procedure. 

Solutions of the Schwinger-Dyson equation are the stationary distributions of 
the complex Langevin equation. We solve equation (129!) following [1]. First 
define 

Zr(j)= j^d<PG{ct>)e^t (30) 
where F is a contour over the complex plane. Inserting this into equation (1291) 
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\ Im O 



TT/6 / 



Re O / 



Fig. 1. The region defined by sin(30) < 0, where 9 is the argument of (p, is shaded. 
Any path starting and ending at infinity within these wedges corresponds to a 
particular solution of zero dimensional icj)"^ theory. 



one gets: 

= - G{<f))e^ 
This equation can be solved for 



G(0) 



-5(0) 



oi0 



(31) 



(32) 



and r is contour that connects the zeros of e ^^'I'^+^'t' on the complex (p plane. 



Now consider polynomial actions, 



i=i ' 



(33) 



The contours will be defined by m wedges, where m is the order of 5(0), 
such that Re{gm4>"^) — +oo as |0| — > oo. Contours obtained by deforming 
r without crossing singularities of e"'^'-''^) and keeping boundary points fixed 
result in the same generating function. Note that set of all Z^lj) will not be 
independent. The Schwinger-Dyson equation will be of order (m — 1), and 



will have (m — 1) independent solutions. For example, if S 



then e-^^"^) 



will have three zeros in the complex plane. Figure [Hwhich will define three 
different generating functions. However the Schwinger-Dyson equation will be 
a second order linear differential equation, which has two independent solu- 
tions. Therefore any two of the three possible paths will define an independent 
solution set. 



Since the Schwinger-Dyson equation is linear, any linear combination of the 
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(independent) solutions will also be a solution, 



Z{j) = J2ar,ZrAj), (34) 

where Tj define an independent subset of paths F. Now defining the distribu- 
tion Pr{(pR) on the real plane as 

= (35) 

which can always be done by a real parametrization of the complex contour 
r, one sees that the equilibrium distribution can be written as a linear com- 
bination 

Pe,(0ij) = E«r,Pr,(0R), (36) 

which is exactly the conjecture that was made by Salcedo [H], where he de- 
rived the same result for a general complex distribution by considering the 
stationary solutions of the pseudo Fokker-Planck equation (l24l) to be real- 
ized as distributions rather than functions. Actually, we managed to refine 
his conjecture (which considers a sum over all F instead of F/ on the RHS of 
equation (1361) ) by showing that not all of Pr{(pR) are independent through the 
use of Schwinger-Dyson equations. A final note is that the coefficients ar may 
depend on initial conditions. We will illustrate these points with numerical 
examples in the next section. |9] has other examples discussed along the lines 
mentioned here. 



3. 1 Zero Dimensional Examples 



We consider two examples here. 



SM = 




In both cases, we will derive an independent set of generating functionals and 
identify the particular solution of the Schwinger-Dyson equation to which the 
simulation converges by observing the sampling points in the complex plane. 
We will see that the change of initial conditions may change the resulting sta- 
tionary distribution. Zero dimensional field theories have been heavily stud- 
ied with complex Langevin equations before, e.g. |26|l3ll27|28|13|29|30ll9] . What 
makes our presentation different from the previous studies is the relation to 
complex path integral solutions of Schwinger-Dyson equations. 

In our simulations, we used the Euler method which is a first-order algorithm, 
see e.g. [16]. The complex Langevin equation (ITTIl with this discretization is 
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Table 1 

Comparison of Langevin simulation results of action Si{(j)) with the correlators of 
the generating functions z[^\j) and Z^\j). The simulation ran from Tj = to 
Tf = 1000 with At = 0.001 with different initial conditions. We start calculating 
expectation values after r = 5. Error bars stand for the standard deviation over 
50 runs. 5 diverging paths discarded for the second initial condition. Average is 
over 45 converging paths. See [31J for a justification of this procedure. The last two 
rows show the corresponding exact values for the generating functions calculated 
by numerical integration. 







m 


= 


-0.0034 - iO.7289 
±0.0190 ± i0.0076 


i0.0053 
±0.0003 ± i0.0364 


0.0025 - fO.9998 
±0.0292 ± i0.0347 


m 


= 5i 


-0.0016 - i0.7315 
±0.0176 ± i0.0079 


i0.0026 
±0.0003 ± i0.0325 


0.0012 - fl.0080 
±0.0260 ± «0.0330 


m 


= l-i 


-0.0049 - iO.7289 
±0.0246 ± i0.0071 


i0.0085 
±0.0003 ± i0.0457 


0.0053 - iO.9994 
±0.0343 ± i0.0307 






-i0.7290 





—i 


7(1) 




-0.6313 ± iO.3645 





—i 



given by 



0i?(rj+i) = (pR{Tj) - Re 

r,+i = r, ± Ar, j E Z, (38) 

where Ar is the time step, and rij is a Gaussian random variable with zero 
mean and unit variance satisfying 

iVj) = 0, iVjVk) = Sjk. (39) 

All our simulations run from r^ = to rj = 1000 with Ar = 0.001. We start 
calculating expectation values after r = 5. Error bars stand for the standard 
deviation of 50 runs. 

For Si{(j)), an independent set of solutions to the Schwinger- Dyson equation 
(1291) can be written by connecting the three zeros shown in Figure [H i.e. too, 
e~*^cxD and e~'^^oo. We choose the following generating functions z[^\j) and 
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dS 
dS 



Ar ± V2ATr]j, 
Ar, 





-0.2 
-0.4 
-0.6 

-e- 

^ -1 
-1.2 
-1.4 
-1.6 



-1.8 



-2.5 




+ 



--Jt +++ + 



+ ++ 



+ 
+ 



++ 



4- 



+ 



-1.5 



-0.5 0.5 
Re (|)(t) 



1.5 



2.5 



Fig. 2. A sample path for the complex Langevin simulation of the action Si ((/>). 
Parameters are the same as in Table [TJ 1000 sample points are shown with equal 
time intervals. Initial conditions is (j){0) = 0. 



/ri,2 d(p exp [-i^ + 
T2 = [e~*^oo,0] + [0,ioo]. 



(40) 



We expect the result of the complex Langevin simulation to converge to a 
linear combination of the distributions defined by these generating functions. 

Table [1] shows the results of simulations for this theory. We compare with the 
exact results for the Green's functions of z[^\j) and Z2^\j). Figure [2] shows 
a sample path for the simulation of this action. We repeated the simulations 
with different initial conditions, some of which are given in Table [H For the 
converging simulations, the sample paths localized around the same region as 
of Figure [2] and the obtained numerical values were similar to those of Table 
[H For some initial conditions with large \dS/d(j)\ values (e.g. 0(0) = 5i), we 
observed nonconverging paths, which could be made to converge (and localized 
in the same region of Figure [2]) by decreasing the step size. This behavior can 
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Table 2 

Comparison of Langevin simulation results of action 5*2 (<?!') with correlators of the 

(2) (2) (2) 

generating functions Za (j), (j) and Z\ (j). The simulation ran form Tj = 
to Tj = 1000 with Ar = 0.001 with different initial conditions. We start calculating 
expectation values after r = 5. Error bars stand for the standard deviation over 
50 runs. The last two rows show the corresponding exact values for the generating 
functions calculated by numerical integration. 





0(0) = 


no convergence 


no convergence 


no convergence 


no convergence 


0(0) = i 


-0.0014 + 10.9781 
±0.0126 ±i0.0057 


-0.6765 - i0.0029 
±0.0078 ± i0.0295 


0.0038 
±0.0452 ±j0.0003 


-1.0014 - i0.0035 
±0.0462 ± i0.0242 


0(0) = -i 


-0.0008 - iO.9794 
±0.0138 ± i0.0071 


-0.6781 + 10.0081 
±0.0098 ± i0.0318 


0.0024 
±0.0480 ±j0.0003 


-1.0054 - i0.0041 
±0.0312 ±j0.0491 


7(2) 


jO.9777 


-0.6760 





-1 


7(2) 


-jO.9777 


-0.6760 





-1 




-0.9777 


0.6760 





-1 



be understood by inspecting the deterministic part of the complex Langevin 
equation (i.e. without a noise term in equation (|T71)). The solution to to the 
deterministic part will be: 

Ur) = j^^, (41) 

1 + l(f)d{0)T 

where 0d(O) is the initial condition. We see that = is a global attractor 
for all points except the positive imaginary axis. Any path starting from the 
positive imaginary axis will go to infinity staying on the imaginary axis (i.e. 
ioo) in finite time. When the noise term is included, which points along the 
real axis, these diverging paths will come out the positive imaginary axis and 
eventually approach the sampling region shown in Figure [2l However, numeri- 
cally these paths may cause a problem. When the step size is not small enough, 
the simulation may go to infinity around these points in finite time. This is 
called the "runaway solution" problem, see [31] for other examples and more 
details. Despite this numerical problem, which can be cured by smaller step 
sizes, the simulations suggest that for the action S'i(0) the complex Langevin 
algorithm always converges. Inspecting Table [T] we see that the distribution 
defined by z[^\j) has correlators within the error range of numerical data. 
The sample path of Figure [2] shows that the simulation does sample around the 
path of z[^\j). It definitely does not sample around the positive imaginary 
axis. Based on these observations, we conjecture that the complex Langevin 
simulation for Si{(j)) = theory always converge to the distribution defined 

by the generating function z[^\j). 

Now we turn to 5*2(0) = —0^/4. Note that this action is real, however e"'^*^'^-' 
does not define a probability distribution on the real line, as the integral 
is divergent. The complex evolution of the associated Langevin equation is 
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Fig. 3. A sample path for the complex Langevin simulation of the action 5*2 ((/>). 
Parameters are the same as in Table [2l 1000 sample points are shown with equal 
time intervals. Initial conditions is (p{0) = i. 

introduced by a choice of complex initial conditions. This procedure, choosing 
complex initial conditions, will turn every real Langevin process to a complex 
Langevin process. The statements we made about complex Langevin equations 
are valid for these cases also. 



An independent set of generating functions for this theory can be written by 
connecting the four zeros of 5*2 (0) on the complex plane, i.e. e^^oo, e~*?cxD, 
e*^cxD and e~*^oo. We choose the following generating functions Zi^\j), 
Zi'\j) andZf (j): 



^(2) . ^ /r.....^'/'exp{^+J0} 
/r...^0exp{^} ' 
r„= [e'^oo,0] + [0,e*5oo], 
Ffe = [e-'^oo,0] + [0,e-'?oo], 

r,= [e-'^oo,0] + [0,e''roo]. (42) 



We expect to see the complex Langevin simulation converge to a linear com- 
bination of the distributions defined by these generating functions. 
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-0.5 



^ -1.5 - 
E 



-2 




+ *++ +4 + + 

+ t 



-2.5 



-1.5 



-0.5 





Re (M-u) 



0.5 



1.5 



Fig. 4. A sample path for the complex Langevin simulation of the action 5*2 ((/>). 
Parameters are the same as in Table [2l 1000 sample points are shown with equal 
time intervals. Initial conditions is (p{0) = —i. 

Table [2] shows the simulation results with different initial conditions compared 
with the correlators of the generating functions. In contrast to the previous 
case, we see that the results of the simulation is initial value dependent. Initial 
points on the real axis do not converge at all. Initial values above and below 
the real axis converge, but the sample points localize in different regions of 
the complex (p plane giving different results, see Figures [3] and IH We can 
understand this behavior again by inspecting the deterministic part of the 
complex Langevin equation. This time the solution will be: 

Mr) = . ^f°' ,^ (43) 

where 0d(O) is the initial condition and the square root function gives the 
principal root. We see that = is a global attractor for every point on the 
complex plane except the real line. The origin will repell any path starting 
on the real line; these paths will diverge in finite time. When the real noise 
term is added, simulations with real initial points will stay on the real line 
and due to the repulsion they will diverge. Initial points on the upper /lower 
half of the complex plane will be attracted by the origin and the simulations 
will sample in the upper /lower half plane as seen in Figures [3] and IH As a 
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result, upper half plane initial points will converge to a different solution of 
the Schwinger-Dyson equation than those of the lower half plane initial points. 
Table [2] shows two simulations starting with 0(0) = i and 0(0) = —i. We see 
that the correlators of Zj^^ are in the error range of the former initial condition, 

(2) 

while the correlators of are in the range of the latter simulation. Based 
on these observations, we conjecture that for 6*2(0), the complex Langevin 
equation will diverge if the initial condition is real, converge to the distribution 

defined by Zj^\j) if the initial condition is on the upper half of the complex 

(2) 

plane or converge to the distribution defined by (j) otherwise. 

One might suspect that either or both of the two stationary states we dis- 
cussed are quasi-stationary states and consequently expect to see the system 
converge to a unique stable stationary state after a long enough simulation 
time. Assuming convergence to a stationary state, we will argue that this is 
not the case, initial points above and below the real line will behave differ- 
ently in the whole range of simulation. Since the path of the complex field 
is continuous (but not differentiable), a complex Langevin process starting 
from above the real line will never end up below the real line. The reverse 
statement holds also. The reason is that any path going from one half plane 
to the other must pass through the real line, and once the path is on the 
real line, it stays on the real line. Both the noise term and the drift term 
points along the real line. Furthermore, the path on the real line will show 
nonconvergent behavior as discussed above. So either, all paths diverge at the 
end, or, assuming convergence, initial point above and below the real line end 
up sampling in different regions of the complex plane, converging to different 
stationary distributions. This initial value dependence is crucial in complex 
Langevin equations. 

Actions 5*1(0) and 5*2(0) were studied in the context of "PT-symmetric quan- 
tum field theories in [SB] , where only one of the path integral solutions to the 
Schwinger-Dyson equation is considered, e.g. jMI- In fact, the authors of pB], 
Bernard and Savage, provided a formal proof that the complex Langevin equa- 
tion for 5*1(0) should always converge to the distribution defined by z[^\j), 
regardless of the initial condition, which we have also observed in our simula- 
tions. The formal proof is based on the methods of [23] and involves defining 
-P(0c, r) of equation (12^ as a projection of P(0_r, 0/, r) to the complex in- 
tegration contour Fi (after necessary partial integrations), as opposed to the 
real line as of equation fl26|l . Then the associated Fokker-Planck Hamiltonian 
will have a real and positive spectrum with one zero eigenvalue. Bernard and 
Savage argued that with slight modifications the same reasoning applies for 
5*2(0) with the integration contour of Fj,. They also noted that in simulations 
of 5*2(0) one has to choose the initial point 0(0) to be in the lower half of 
the complex plane or else the numerical simulations will be unstable. In 
our studies, we observed instabilities for initial points on the real line, which 
we interpreted to be nonconvergent behavior. Initial points in the lower and 
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upper halves of the complex plane were observed to converge to different 
probability distributions (solutions of the Schwinger-Dyson equation). Our re- 
sults suggest that the proof given by Bernard and Savage should be modified 
to include the effects of different initial conditions and all possible solutions of 
the Schwinger-Dyson equation. We believe that more attention must be paid 
to the boundary conditions of P{(f>R, (pi, r) when doing the projection. 

Another point to note is that in both cases, we could not recover the whole 
solution set of the Schwinger-Dyson equations from the complex Langevin 
equation. In general, 0^ theories will have — 1 independent solutions to the 
Schwinger-Dyson equations. In the cases that we studied, we could recover 
only — 2 of them. We may need to consider other stochastic processes to 
recover the whole solution set. 



4 Lattice Study 



It is trivial to generalize this discussion to a lattice. In particular, consider a 
general Euclidean scalar field theory with a polynomial potential term, 

C = d^(f){x) + J2j<f)'{x). (44) 
1=1 '■ 

For simplicity we consider a lattice in one dimension, which we call time and 
denote by t. The generating functional of the theory is now a function of m 
variables, where m is the number of lattice points. We denote the generating 
functional by Z{j) = Z{ji, . . . ,jm)- ji is the source at i*^ lattice point. We 
will use j = to mean ji = . . . = jm = 0. By definition, correlation functions 
are given by 



(45) 



Next, we introduce the difference operators on the lattice, 

A+ = (f)i+i - (f)i_i, A_ = (pi- (pi^i, (46) 

where (pi stands for the field value at i^^ lattice point. The Schwinger-Dyson 
equation for the scalar field theory fj44|) can be written by a coupled set of 
partial differential equations. Using a centered discretization for time deriva- 
tives, the system is composed of a partial differential equation for every lattice 
point, given by 

(a.a_|+|:«,|;^)zw=..zo). (47) 



17 



where the lattice spacing is set to one. 



This discretization requires one to set boundary conditions on Z{j) at initial 
and final lattice points. There are many possible choices, we choose periodic 
boundary conditions. We assume that the whole space is filled with a lattice 
of period m. 

In the absence of the kinetic term, each lattice point acts independently. There- 
fore, study of the complex Langevin equation for zero dimensional case im- 
mediately tells us what will happen on the lattice. When the kinetic term is 
included, couplings between lattice points take action and problem is more 
complicated. 

A general solution to the lattice Schwinger-Dyson equation can be written as 
^(^) = T7 / f d(j)m 

- E i^-^^^^+^-<P^ + E J^'j + E ' (48) 

where Fj denote contours that connect the zeros of the integrand on complex 
(pi plane. We again normalize so that Z{0) = 1. For each lattice point there are 
n independent solutions, which leads to m x n independent solutions to the 
whole lattice Schwinger-Dyson equation. Because of linearity any combination 
of these solutions is also a solution. We note that there is no ambiguity in the 
solution written in this form, one can do the integrations in any order. 

We turn to the complex Langevin equation for this problem. We again in- 
troduce a fictitious time coordinate r. The complex Langevin system is now 
written in terms of stochastic variables (pjir) = (pRjij) + icpijij). For each 
lattice point j = 1, . . . ,m, there is a stochastic equation: 

n 

d<P,{r) = -A+A_0,(r) -Y.giU'^)' + dw,{T), (49) 

1=1 

where dwj{T) are m independent Wiener processes normalized as before. A 
repetition of the analysis of the previous section is sufficient to conclude that 
the stationary distributions of this set of stochastic equations will satisfy lat- 
tice Schwinger-Dyson equations (j4j 



Let's consider again the —0^/4 theory, this time in one dimension. The La- 
grangian is given by 

We again note that this theory is bottomless, a normal (real line contour) 
path integral solution to the Schwinger-Dyson equations does not exist. How- 
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Table 3 

A complex Langevin simulation is done for the one dimensional —0^/4 theory. 128 
lattice points are taken with spacing set to 1. An Euler method is used for fictitous 
time evolution. Simulations ran from Tj = to = 10000, with Ar = 0.001. 
The table below lists the first four equal time correlation functions obtained for 
two different initial conditions. We observed translation symmetry on the equal 
time correlation functions as expected. We list average values over lattice points, 
e.g. ((/>) = Si=i {<Pi)- Standard deviations are calculated over lattice points. We 
also considered other initial conditions, which resulted in either nonconvergence 
or convergence to one of the two cases shown below. A detailed study of initial 
condition dependence is beyond the scope of this paper. 





W 


if') 


(03) 


(<^0 


</>(0,t) = -i 


-j0.8891 


-0.5537 + i0.0001 


0.0001 


-0.6307 - iO.OOOl 




±0.0032 ±i0.0012 


±0.0013 ±j0.0065 


±0.0083 ± i0.0018 


±0.0069 ± j0.0063 


Ho,t) = i 


j0.8891 
±0.0037 ±i0.0013 


-0.5537 - iO.OOOl 
±0.0013 ± j0.0075 


0.0001 
±0.0095 ± i0.0021 


-0.6306 + iO.OOOl 
±0.0071 ± j0.0071 



ever, complex contour contour solutions do exist. Table [3] shows the results 
of numerical simulations for this theory on a one dimensional lattice. As in 
zero dimensional case, we see that the complex Langevin equation has at least 
two different stationary distributions. Choice of initial conditions can alter the 
stationary distributions. We note that this theory was also studied in [26] and 
there initial conditions were restricted to lower half of the complex plane. This 
led to the observation of only one of the stationary distributions, which was 
concluded to be the stationary distribution that led to a PT-symmetric theory. 
This could be described by choosing Tb of the previous section at each lattice 
point as the integration contour for the generating functional (HHl) . Table [3] 
suggests that the other solution is given by choosing at every point. 

More studies on the lattice must be done to understand the convergence behav- 
ior of complex Langevin equations. Here, as well as initial conditions, boundary 
conditions on the lattice may also take effect. Some other specific questions are 
Usted in the next section. We leave detailed lattice studies to a coming paper. 
The point of this section is to demonstrate that the results of the previous 
section applies to the lattice as well. 



5 Discussion and Conclusion 



Salcedo [9] suggested that stationary distributions of the complex Langevin 
equation may be interpreted as different phases of the associated quantum 
field theory. Some authors used Langevin and complex Langevin equation 
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to give physical meaning to bottomless actions, e.g. [33|28f34ll35fl3j^l . Our 
discussion shows that a rephrasal of these questions is to ask which solutions 
of the Schwinger-Dyson equations define a phase of quantum field theory. Here 
we discuss this point, and the relation of different phases to the boundary 
conditions of Schwinger-Dyson equations |T|2]. 

Schwinger-Dyson equations are differential equations and admit more than 
one solution. Therefore it is necessary to set boundary conditions to specify 
the particular solution one is looking for. If one is solving for a quantum field 
theory using Schwinger-Dyson equations, it seems reasonable to choose the 
boundary condition so that the solution is the standard path integral over real 
fields. However, in many cases this solution actually will not be the physical 
one, e.g. symmetry breaking phases. Also, in theories with actions unbounded 
below, the integrals over real fields are not even convergent. In these cases, it 
is reasonable to look at other solutions of the Schwinger-Dyson equations and 
study them as possible generating functionals of the associated quantum field 
theory. Of course, different solutions require specification of different bound- 
ary conditions. One way of specifying different solutions is to consider path 
integrals over complex paths (as opposed to real paths) that connect zeros of 
the partition function on the complex plane, as was demonstrated in equa- 
tion (IHT]) . Some of these also happen to be the stationary distributions of the 
associated complex Langevin equation constructed by Salcedo |9] as shown 
in this paper. The problem with this approach is the large number of differ- 
ent boundary conditions/solutions and if all these different solutions define a 
phase or vacuum of the associated quantum field theory. A possible reduction 
of the solution set comes from taking the thermodynamic limit of the lattice. 
These issues are discussed in detail in [TP] . 

In the light of the discussion above, one concludes that to study the phase 
structure of quantum field theories, one needs to study the different solutions 
to Schwinger-Dyson equations. This task requires new numerical methods in 
the study of quantum field theory. One suggestion is the Source- Galerkin 
method, see e.g. [3B] and references therein. This method proposes an ex- 
pansion of the generating functional in polynomials of the source term and 
optimizes this expansion by a Galerkin procedure using the Schwinger-Dyson 
equation. It is successful in many problems, but also proved to be very difficult 
in many other cases. Another approach is given by mollification of the path in- 
tegral weight [SZj. The connection between complex Langevin and Schwinger- 
Dyson equations suggests the use of complex Langevin simulations as a nu- 
merical method to study the phase structure of quantum field theories. 



Some of these references use nonconstant kernels in the complex Langevin equa- 
tion which enlarges the set of stationary distributions [9J and changes the Schwinger- 
Dyson equation, e.g. |13j 
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Some questions and speculations in this quest are: 

• Is it possible to know a priori if the complex Langevin equation will con- 
verge? 

• What is the exact relation between the initial condition and the stationary 
distributions of the complex Langevin equation? 

• Is it possible to recover the whole solution set of the Schwinger-Dyson 
equation using the complex Langevin equation? If not, can we use other 
stochastic systems to recover the whole set? Also, what is special about the 
recovered solutions? 

• How are these results modified in the continuum? There are an infinite 
number of solutions to Schwinger-Dyson equations on the lattice. The con- 
tinuum limit may cause a collapse in the solution set. This is definitely the 
case if space time translational invariance is required. [II2]- Can one see this 
collapse using complex Langevin equations? This will be very important in 
understanding the phase structure of quantum field theories. 

We address some of these problems in future work. 
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